COVID19 GER

Prepare data

Read and format data

df_ger_covid <- read_csv('Germany_timeseries_prep.csv')
df_ger_socdist <- read_csv('Germany_socdist_fb_kreis.csv')
df_ger_ctrl <- read_delim('Germany_controls.csv', delim = ';')

# prevalence 
df_ger_covid_clean <- df_ger_covid %>% mutate(date = as.Date(date, "%d%b%Y"),
                                  kreis = as.character(kreis)) %>% 
  filter(date >= '2020-03-09' & date <= '2020-03-31') %>% 
  group_by(kreis) %>% 
  mutate(time = row_number()) %>% 
  ungroup() %>% 
  dplyr::select(-runday, -kreis_name, -ewz, -shape__area, 
                -cumcase, -anzahlfall, -popdens) %>%
  dplyr::rename(pers_o = open, 
         pers_c = sci,
         pers_e = extra,
         pers_a = agree,
         pers_n = neuro)

df_ger_covid_clean %>% head()

# social distancing
df_ger_socdist_clean <- df_ger_socdist %>% 
  filter(date >= '2020-03-09' & date <= '2020-03-31') %>% 
  mutate(kreis = as.character(kreis)) %>%
  group_by(kreis) %>% 
  mutate(time = row_number()) %>% 
  ungroup() %>% 
  rename(socdist_single_tile = all_day_ratio_single_tile_users) %>%
  select(kreis, time, socdist_single_tile)

df_ger_socdist_clean %>% head()

NA

# controls 
df_ger_ctrl_clean <- df_ger_ctrl %>% select(-kreis_nme) %>%
    mutate(kreis = as.character(kreis),
           hospital_beds = as.numeric(str_replace(hospital_beds, ',', '.')))

df_ger_ctrl_clean %>% head()

NA
NA

# merge
df_ger <- df_ger_covid_clean %>% 
  plyr::join(df_ger_socdist_clean, by = c('kreis', 'time')) %>% 
  inner_join(df_ger_ctrl_clean, by = 'kreis')

df_ger

NA

Explore data

Plot prevalence over time

df_ger %>% ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=kreis, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall prevalence over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_ger %>% mutate(prev_tail = cut(.[[i]], 
      breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
      labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(prev_tail != 'center') %>%
  ggplot(aes(x=time, y=rate_day)) + 
  geom_point(aes(col=kreis, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~prev_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

Plot social distancing (single tile) over time

df_ger %>% ggplot(aes(x=time, y=socdist_single_tile)) + 
  geom_point(aes(col=kreis, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall social distancing (single tile) over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_ger %>% mutate(prev_tail = cut(.[[i]], 
      breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
      labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(prev_tail != 'center') %>%
  ggplot(aes(x=time, y=socdist_single_tile)) + 
  geom_point(aes(col=kreis, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~prev_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

Control for weekend effect

weekend <- c(6, 7, 13, 14, 20, 21)

df_ger_loess <- df_ger %>% filter(!time %in% weekend) %>% 
  split(.$kreis) %>%
  map(~ loess(socdist_single_tile ~ time, data = .)) %>%
  map(predict, 1:23) %>% 
  bind_rows() %>% 
  gather(key = 'kreis', value = 'loess') %>% 
  group_by(kreis) %>% 
  mutate(time = row_number())


df_ger <- df_ger %>% merge(df_ger_loess, by=c('kreis', 'time')) %>% 
  mutate(socdist_single_tile_clean = ifelse(time %in% weekend, loess,
                                            socdist_single_tile)) %>%
  arrange(kreis, time)


df_ger %>% ggplot(aes(x=time, y=loess, group=kreis)) +
  geom_line()

NA
NA

df_ger %>% ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=kreis, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  theme(legend.position="none") +
  ggtitle("Overall social distancing (single tile) over time")

pers <- c('pers_o', 'pers_c', 'pers_e', 'pers_a', 'pers_n')

for (i in pers){

gg <- df_ger %>% mutate(prev_tail = cut(.[[i]], 
      breaks = c(-Inf, quantile(.[[i]], 0.2), quantile(.[[i]], 0.8), Inf),
      labels = c('lower tail', 'center', 'upper tail'))) %>% 
  filter(prev_tail != 'center') %>%
  ggplot(aes(x=time, y=socdist_single_tile_clean)) + 
  geom_point(aes(col=kreis, size=popdens)) + 
  geom_smooth(method="loess", se=T) + 
  facet_wrap(~prev_tail) + 
  theme(legend.position="none") +
  ggtitle(i)

print(gg)
}

Correlations

df_ger %>% group_by(kreis) %>% 
  summarize_if(is.numeric, mean) %>% 
  select(-kreis, -time) %>% 
  cor(use='pairwise.complete') %>% 
  round(3)

                          pers_e pers_a pers_c pers_n pers_o rate_day socdist_single_tile  women academics    cdu    afd hospital_beds
pers_e                     1.000  0.223  0.255 -0.375  0.277    0.168               0.149 -0.011     0.162  0.111 -0.108        -0.037
pers_a                     0.223  1.000  0.347 -0.379  0.167    0.135               0.002 -0.009     0.180  0.023  0.015         0.004
pers_c                     0.255  0.347  1.000 -0.373 -0.063   -0.010               0.030 -0.046    -0.045 -0.006  0.140        -0.157
pers_n                    -0.375 -0.379 -0.373  1.000 -0.046   -0.127              -0.141 -0.012    -0.068 -0.061  0.064         0.111
pers_o                     0.277  0.167 -0.063 -0.046  1.000    0.115               0.090 -0.094     0.431 -0.053 -0.187         0.281
rate_day                   0.168  0.135 -0.010 -0.127  0.115    1.000               0.265 -0.026     0.128  0.136 -0.135        -0.060
socdist_single_tile        0.149  0.002  0.030 -0.141  0.090    0.265               1.000  0.049    -0.018  0.149 -0.301        -0.243
women                     -0.011 -0.009 -0.046 -0.012 -0.094   -0.026               0.049  1.000    -0.001  0.074 -0.089         0.009
academics                  0.162  0.180 -0.045 -0.068  0.431    0.128              -0.018 -0.001     1.000 -0.208 -0.117         0.335
cdu                        0.111  0.023 -0.006 -0.061 -0.053    0.136               0.149  0.074    -0.208  1.000 -0.157        -0.091
afd                       -0.108  0.015  0.140  0.064 -0.187   -0.135              -0.301 -0.089    -0.117 -0.157  1.000        -0.030
hospital_beds             -0.037  0.004 -0.157  0.111  0.281   -0.060              -0.243  0.009     0.335 -0.091 -0.030         1.000
tourism_beds              -0.115 -0.134 -0.102  0.007 -0.108   -0.038              -0.012  0.003    -0.156  0.115 -0.011        -0.043
gdp                        0.121  0.064  0.000 -0.067  0.247    0.110              -0.047 -0.044     0.374 -0.078 -0.133         0.266
manufact                  -0.009 -0.061 -0.051 -0.007 -0.038    0.076              -0.071  0.026    -0.088  0.019  0.065         0.070
airport                   -0.188 -0.167 -0.112  0.171 -0.222   -0.031              -0.142  0.043    -0.331  0.105  0.183        -0.046
age                        0.007  0.086  0.008 -0.016 -0.063   -0.034              -0.017 -0.046    -0.051  0.020  0.091         0.029
popdens                    0.062  0.038  0.018 -0.080  0.072    0.092               0.110 -0.001     0.056 -0.016 -0.028         0.024
loess                      0.139 -0.013  0.020 -0.143  0.105    0.229               0.979  0.040     0.012  0.124 -0.344        -0.229
socdist_single_tile_clean  0.138 -0.015  0.020 -0.142  0.104    0.228               0.979  0.041     0.010  0.123 -0.344        -0.228
                          tourism_beds    gdp manufact airport    age popdens  loess socdist_single_tile_clean
pers_e                          -0.115  0.121   -0.009  -0.188  0.007   0.062  0.139                     0.138
pers_a                          -0.134  0.064   -0.061  -0.167  0.086   0.038 -0.013                    -0.015
pers_c                          -0.102  0.000   -0.051  -0.112  0.008   0.018  0.020                     0.020
pers_n                           0.007 -0.067   -0.007   0.171 -0.016  -0.080 -0.143                    -0.142
pers_o                          -0.108  0.247   -0.038  -0.222 -0.063   0.072  0.105                     0.104
rate_day                        -0.038  0.110    0.076  -0.031 -0.034   0.092  0.229                     0.228
socdist_single_tile             -0.012 -0.047   -0.071  -0.142 -0.017   0.110  0.979                     0.979
women                            0.003 -0.044    0.026   0.043 -0.046  -0.001  0.040                     0.041
academics                       -0.156  0.374   -0.088  -0.331 -0.051   0.056  0.012                     0.010
cdu                              0.115 -0.078    0.019   0.105  0.020  -0.016  0.124                     0.123
afd                             -0.011 -0.133    0.065   0.183  0.091  -0.028 -0.344                    -0.344
hospital_beds                   -0.043  0.266    0.070  -0.046  0.029   0.024 -0.229                    -0.228
tourism_beds                     1.000 -0.094   -0.051   0.297  0.075  -0.157  0.031                     0.031
gdp                             -0.094  1.000    0.339  -0.085 -0.014   0.082 -0.026                    -0.027
manufact                        -0.051  0.339    1.000   0.136 -0.049   0.072 -0.119                    -0.120
airport                          0.297 -0.085    0.136   1.000  0.077  -0.179 -0.169                    -0.169
age                              0.075 -0.014   -0.049   0.077  1.000  -0.038 -0.021                    -0.021
popdens                         -0.157  0.082    0.072  -0.179 -0.038   1.000  0.094                     0.094
loess                            0.031 -0.026   -0.119  -0.169 -0.021   0.094  1.000                     1.000
socdist_single_tile_clean        0.031 -0.027   -0.120  -0.169 -0.021   0.094  1.000                     1.000

Modelling

Prepare functions

# function calculates all relevant models
run_models <- function(y, lvl1_x, lvl2_x, lvl2_id, data, ctrls=F){

  # subset data
  data = data %>% 
    dplyr::select(all_of(y), all_of(lvl1_x), all_of(lvl2_x), all_of(lvl2_id), 
                  popdens, rate_day)
  data = data %>% 
    dplyr::rename(y = all_of(y),
           lvl1_x = all_of(lvl1_x),
           lvl2_x = all_of(lvl2_x),
           lvl2_id = all_of(lvl2_id)
           )
  
  # configure optimization procedure
  ctrl_config <- lmeControl(opt = 'optim', maxIter = 100, msMaxIter = 100)

  # baseline
  baseline <- lme(fixed = y ~ 1, random = ~ 1 | lvl2_id, 
                    data = data,
                    correlation = corAR1(),
                  control = ctrl_config,
                  method = 'ML')

  # random intercept fixed slope
  random_intercept <- lme(fixed = y ~ lvl1_x + lvl2_x, 
                          random = ~ 1 | lvl2_id,
                            data = data,
                            correlation = corAR1(),
                  control = ctrl_config,
                  method = 'ML')

  # random intercept random slope
  random_slope <- lme(fixed = y ~ lvl1_x + lvl2_x, 
                      random = ~ lvl1_x | lvl2_id, 
                        data = data,
                        correlation = corAR1(),
                  control = ctrl_config,
                  method = 'ML')

  # cross level interaction
  interaction <- lme(fixed = y ~ lvl1_x * lvl2_x, 
                     random = ~ lvl1_x | lvl2_id, 
                       data = data,
                       correlation = corAR1(),
                  control = ctrl_config,
                  method = 'ML')
  
  # create list with results
  results <- list('baseline' = baseline, 
                  "random_intercept" = random_intercept, 
                  "random_slope" = random_slope,
                  "interaction" = interaction)
  
  
  if (ctrls == 'dem' | ctrls == 'prev'){
    
    # random intercept random slope
    random_slope_ctrl_dem <- lme(fixed = y ~ lvl1_x + lvl2_x + popdens,
                              random = ~ lvl1_x | lvl2_id, 
                          data = data,
                          correlation = corAR1(),
                    control = ctrl_config,
                  method = 'ML')
  
    # cross level interaction
    interaction_ctrl_main_dem <- lme(fixed = y ~ lvl1_x * lvl2_x + popdens,
                             random = ~ lvl1_x | lvl2_id, 
                         data = data,
                         correlation = corAR1(),
                    control = ctrl_config,
                  method = 'ML')
  
    # cross level interaction
    interaction_ctrl_int_dem <- lme(fixed = y ~ lvl1_x * lvl2_x + lvl1_x * popdens,
                             random = ~ lvl1_x | lvl2_id, 
                         data = data,
                         correlation = corAR1(),
                    control = ctrl_config,
                  method = 'ML')        
    
    # create list with results
    results <- list('baseline' = baseline, 
                    "random_intercept" = random_intercept, 
                    "random_slope" = random_slope,
                    "interaction" = interaction,
                    "random_slope_ctrl_dem" = random_slope_ctrl_dem,
                    "interaction_ctrl_main_dem" = interaction_ctrl_main_dem,
                    "interaction_ctrl_int_dem" = interaction_ctrl_int_dem)
  }
  
  if (ctrls == 'prev'){
  
    # random intercept random slope
    random_slope_ctrl_prev <- lme(fixed = y ~ lvl1_x + lvl2_x + popdens + rate_day,
                              random = ~ lvl1_x + rate_day | lvl2_id, 
                          data = data,
                          correlation = corAR1(),
                          control = ctrl_config,
                  method = 'ML')  
    
        # cross level interaction
    interaction_ctrl_main_prev <- lme(fixed = y ~ lvl1_x * lvl2_x + popdens + rate_day,
                             random = ~ lvl1_x | lvl2_id, 
                         data = data,
                         correlation = corAR1(),
                    control = ctrl_config,
                  method = 'ML')
  
  
    # cross level interaction
    interaction_ctrl_int_prev<- lme(fixed = y ~ lvl1_x * lvl2_x + lvl1_x * popdens + rate_day,
                             random = ~ lvl1_x + rate_day | lvl2_id, 
                         data = data,
                         correlation = corAR1(),
                          control = ctrl_config,
                  method = 'ML')
  
    # create list with results
    results <- list('baseline' = baseline, 
                    "random_intercept" = random_intercept, 
                    "random_slope" = random_slope,
                    "interaction" = interaction,
                    "random_slope_ctrl_dem" = random_slope_ctrl_dem,
                    "interaction_ctrl_main_dem" = interaction_ctrl_main_dem,
                    "interaction_ctrl_int_dem" = interaction_ctrl_int_dem,                    
                    "random_slope_ctrl_prev" = random_slope_ctrl_prev,
                    "interaction_ctrl_main_prev" = interaction_ctrl_main_prev,
                    "interaction_ctrl_int_prev" = interaction_ctrl_int_prev)
  }
  
  if(ctrls == 'exp'){
    # random intercept random slope
  random_slope_exp <- lme(fixed = y ~ (lvl1_x + I(lvl1_x^2)) + lvl2_x, 
                      random = ~ (lvl1_x + I(lvl1_x^2)) | lvl2_id, 
                        data = data,
                        correlation = corAR1(),
                  control = ctrl_config,
                  method = 'ML')

  # cross level interaction
  interaction_exp <- lme(fixed = y ~ (lvl1_x + I(lvl1_x^2)) * lvl2_x, 
                     random = ~ (lvl1_x + I(lvl1_x^2)) | lvl2_id, 
                       data = data,
                       correlation = corAR1(),
                  control = ctrl_config,
                  method = 'ML')  
  
  
  # create list with results
  results <- list('baseline' = baseline, 
                  "random_intercept" = random_intercept, 
                  "random_slope" = random_slope,
                  "interaction" = interaction,                  
                  "random_slope_exp" = random_slope_exp,
                  "interaction_exp" = interaction_exp)
  }
  
  return(results)
        
}

# extracts table with coefficients and tests statistics
extract_results <- function(models) {
  
  models_summary <- models %>% 
  map(summary) %>% 
  map("tTable") %>% 
  map(as.data.frame) %>% 
  map(round, 10) 
  # %>% map(~ .[str_detect(rownames(.), 'Inter|lvl'),])
  
  return(models_summary)
  
}


compare_models <- function(models) {

  mdl_names <- models %>% names()
  
  str = ''
  for (i in mdl_names){
    
    mdl_str <- paste('models$', i, sep = '')
    
    if(i == 'baseline'){
      str <- mdl_str
    }else{
    str <- paste(str, mdl_str, sep=', ')
    }
  }
  
  anova_str <- paste0('anova(', str, ')')
  mdl_comp <- eval(parse(text=anova_str))
  rownames(mdl_comp) = mdl_names
  return(mdl_comp)
}

Rescale Data

df_ger_scaled <- df_ger %>% dplyr::select(kreis, time, pers_o, 
                                 pers_c, pers_e, pers_e, pers_a, pers_n,
                                 popdens, rate_day, socdist_single_tile) %>%
  mutate_at(vars(-kreis, -time), scale)

df_ger_scaled %>% head()

Predict prevalence

prevalence ~ openness

models_o_covid <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_o', 
                         lvl2_id = 'kreis', 
                         data = df_ger_scaled,
                         ctrls = 'dem')

extract_results(models_o_covid)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

compare_models(models_o_covid)

NA

prevalence ~ conscientiousness

models_c_covid <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_c', 
                         lvl2_id = 'kreis', 
                         data = df_ger_scaled,
                         ctrls = 'dem')

extract_results(models_c_covid)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

compare_models(models_c_covid)

NA
NA

prevalence ~ extraversion

models_e_covid <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_e', 
                         lvl2_id = 'kreis', 
                         data = df_ger_scaled,
                         ctrls = 'dem')

extract_results(models_e_covid)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

compare_models(models_e_covid)

NA
NA

prevalence ~ agreeableness

models_a_covid <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_a', 
                         lvl2_id = 'kreis', 
                         data = df_ger_scaled,
                         ctrls = 'dem')

extract_results(models_a_covid)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

compare_models(models_a_covid)

NA
NA

prevalence ~ neuroticism

models_n_covid <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_n', 
                         lvl2_id = 'kreis', 
                         data = df_ger_scaled,
                         ctrls = 'dem')

extract_results(models_n_covid)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

compare_models(models_n_covid)

NA
NA

Predict social distancing

social distancing ~ openness

models_o_socdist <-run_models(y = 'socdist_single_tile', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_o', 
                         lvl2_id = 'kreis', 
                         data = df_ger_scaled,
                         ctrls = 'dem')

extract_results(models_o_socdist)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

compare_models(models_o_socdist)

NA

social distancing ~ conscientiousness

models_c_socdist <-run_models(y = 'socdist_single_tile', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_c', 
                         lvl2_id = 'kreis', 
                         data = df_ger_scaled,
                         ctrls = 'dem')

extract_results(models_c_socdist)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

compare_models(models_c_socdist)

NA
NA
NA

social distancing ~ extraversion

models_e_socdist <-run_models(y = 'socdist_single_tile', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_e', 
                         lvl2_id = 'kreis', 
                         data = df_ger_scaled,
                         ctrls = 'dem')

extract_results(models_e_socdist)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

compare_models(models_e_socdist)

NA
NA
NA

social distancing ~ agreeableness

models_a_socdist <-run_models(y = 'socdist_single_tile', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_a', 
                         lvl2_id = 'kreis', 
                         data = df_ger_scaled,
                         ctrls = 'dem')

extract_results(models_a_socdist)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

compare_models(models_a_socdist)

NA
NA
NA

social distancing ~ neuroticism

models_n_socdist <-run_models(y = 'socdist_single_tile', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_n', 
                         lvl2_id = 'kreis', 
                         data = df_ger_scaled,
                         ctrls = 'dem')

extract_results(models_n_socdist)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_ctrl_dem

$interaction_ctrl_main_dem

$interaction_ctrl_int_dem

compare_models(models_n_socdist)

NA
NA

Explore quadratic trends

prevalence ~ openness

models_o_covid_exp <-run_models(y = 'rate_day',
                         lvl1_x = 'time',
                         lvl2_x = 'pers_o',
                         lvl2_id = 'kreis',
                         data = df_ger_scaled,
                         ctrls = 'exp')

extract_results(models_o_covid_exp)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_exp

$interaction_exp

compare_models(models_o_covid_exp)

NA

prevalence ~ conscientiousness

models_c_covid_exp <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_c', 
                         lvl2_id = 'kreis', 
                         data = df_ger_scaled,
                         ctrls = 'exp')

extract_results(models_c_covid_exp)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_exp

$interaction_exp

compare_models(models_c_covid_exp)

NA

prevalence ~ extraversion

models_e_covid_exp <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_e', 
                         lvl2_id = 'kreis', 
                         data = df_ger_scaled,
                         ctrls = 'exp')

extract_results(models_e_covid_exp)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_exp

$interaction_exp

compare_models(models_e_covid_exp)

NA

prevalence ~ agreeableness

models_a_covid_exp <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_a', 
                         lvl2_id = 'kreis', 
                         data = df_ger_scaled,
                         ctrls = 'exp')

extract_results(models_a_covid_exp)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_exp

$interaction_exp

compare_models(models_a_covid_exp)

NA

prevalence ~ neuroticism

models_n_covid_exp <-run_models(y = 'rate_day', 
                         lvl1_x = 'time', 
                         lvl2_x = 'pers_n', 
                         lvl2_id = 'kreis', 
                         data = df_ger_scaled,
                         ctrls = 'exp')

extract_results(models_n_covid_exp)

$baseline

$random_intercept

$random_slope

$interaction

$random_slope_exp

$interaction_exp

compare_models(models_n_covid_exp)

NA

Create overview table

Define function to create overview tables

summary_table <- function(models, dv_name){

  temp_df_ctrl_main <- NULL
  temp_df_ctrl_int <- NULL
  
  for (i in models){
    results <- i %>% extract_results()
    
    results_ctrl_main <- results$interaction_ctrl_main_dem['lvl1_x:lvl2_x',]
    temp_df_ctrl_main <- temp_df_ctrl_main %>% rbind(results_ctrl_main)
    
    results_ctrl_int <- results$interaction_ctrl_int_dem['lvl1_x:lvl2_x',]
    temp_df_ctrl_int <- temp_df_ctrl_int %>% rbind(results_ctrl_int)
  }
  
  names_ctrl_main <- paste0(dv_name, '~', c('o', 'c', 'e', 'a', 'n'), '*time', '_crtl_popdens')
  names_ctrl_int <- paste0(dv_name, '~', c('o', 'c', 'e', 'a', 'n'), '*time', '_crtl_popdens*time')

  rownames(temp_df_ctrl_main) <- names_ctrl_main
  rownames(temp_df_ctrl_int) <- names_ctrl_int
  
  sum_tab <- rbind(temp_df_ctrl_main, temp_df_ctrl_int) %>% round(4)
  
  return(sum_tab)

} 

Create overview tables

# prevalence
models_prev <- list(models_o_covid, 
                    models_c_covid, 
                    models_e_covid, 
                    models_a_covid, 
                    models_n_covid)

sum_tab_prev <- summary_table(models_prev, dv_name = 'prev')

write.table(sum_tab_prev, quote=F)

Value Std.Error DF t-value p-value
prev~o*time_crtl_popdens 0.0071 0.0029 8798 2.4364 0.0149
prev~c*time_crtl_popdens 3e-04 0.003 8798 0.0944 0.9248
prev~e*time_crtl_popdens 0.0099 0.0029 8798 3.402 7e-04
prev~a*time_crtl_popdens 0.0091 0.0029 8798 3.1314 0.0017
prev~n*time_crtl_popdens -0.0072 0.0029 8798 -2.4642 0.0138
prev~o*time_crtl_popdens*time 0.0069 0.0029 8797 2.3439 0.0191
prev~c*time_crtl_popdens*time 2e-04 0.0029 8797 0.0687 0.9452
prev~e*time_crtl_popdens*time 0.0097 0.0029 8797 3.3233 9e-04
prev~a*time_crtl_popdens*time 0.009 0.0029 8797 3.086 0.002
prev~n*time_crtl_popdens*time -0.0069 0.0029 8797 -2.3608 0.0183

# social distancing
models_socdist <- list(models_o_socdist, 
                       models_c_socdist, 
                       models_e_socdist, 
                       models_a_socdist, 
                       models_n_socdist)

sum_tab_socdist <- summary_table(models_socdist, dv_name = 'socdist')

write.table(sum_tab_socdist, quote=F)

Value Std.Error DF t-value p-value
socdist~o*time_crtl_popdens 0.0041 0.0014 8798 3.0071 0.0026
socdist~c*time_crtl_popdens 0 0.0014 8798 -0.0081 0.9935
socdist~e*time_crtl_popdens 0.003 0.0014 8798 2.1539 0.0313
socdist~a*time_crtl_popdens 0.0017 0.0014 8798 1.2497 0.2114
socdist~n*time_crtl_popdens -0.0015 0.0014 8798 -1.0885 0.2764
socdist~o*time_crtl_popdens*time 0.0041 0.0014 8797 2.9445 0.0032
socdist~c*time_crtl_popdens*time 0 0.0014 8797 -0.0257 0.9795
socdist~e*time_crtl_popdens*time 0.0029 0.0014 8797 2.0978 0.036
socdist~a*time_crtl_popdens*time 0.0017 0.0014 8797 1.2133 0.2251
socdist~n*time_crtl_popdens*time -0.0014 0.0014 8797 -1.0139 0.3106

Conditional random forest analysis

Extract slopes prevalence

# slope prevalence
df_ger_slope_prev <- df_ger %>% split(.$kreis) %>% 
  map(~ lm(rate_day ~ time, data = .)) %>%
  map(coef) %>% 
  map_dbl('time') %>% 
  as.data.frame() %>% 
  rownames_to_column('kreis') %>% 
  rename(slope_prev = '.')

# merge with control variables 
df_ger_slope_prev <- df_ger %>% select(-time, -date, -rate_day, -socdist_single_tile) %>%
  distinct() %>% 
  inner_join(df_ger_slope_prev, by = 'kreis') %>%
  drop_na()

df_ger_slope_prev

NA

Extract slopes social distancing

# slope socdist
df_ger_slope_socdist <- df_ger %>% split(.$kreis) %>%
  map(~ lm(socdist_single_tile ~ time, data = .)) %>%
  map(coef) %>%
  map_dbl('time') %>%
  as.data.frame() %>%
  rownames_to_column('kreis') %>%
  rename(slope_socdist = '.')


# merge with control variables 
df_ger_slope_socdist <- df_ger %>% 
  select(-time, -date, -socdist_single_tile, -rate_day) %>%
  distinct() %>%
  inner_join(df_ger_slope_socdist, by = 'kreis') %>%
  drop_na()

df_ger_slope_socdist

NA

Explore distribution of slopes

df_ger_slope_prev %>% ggplot(aes(slope_prev)) + geom_histogram(bins = 100)

df_ger_slope_socdist %>% ggplot(aes(slope_socdist)) + geom_histogram(bins = 100)

NA
NA

CRF prevalence ~ openness

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_o_fit_prev <- cforest(slope_prev ~ pers_o + women + academics +
                          cdu + afd + hospital_beds + gdp + manufact +
                          airport + age + popdens, 
                         df_ger_slope_prev[-1], 
                         controls = ctrls)

crf_o_varimp_prev <- varimp(crf_o_fit_prev, nperm = 1)
crf_o_varimp_cond_prev <- varimp(crf_o_fit_prev, conditional = T, nperm = 1)

CRF prevalence ~ conscientiousness

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_c_fit_prev <- cforest(slope_prev ~ pers_c + women + academics +
                          cdu + afd + hospital_beds + gdp + manufact +
                          airport + age + popdens, 
                         df_ger_slope_prev[-1], 
                         controls = ctrls)

crf_c_varimp_prev <- varimp(crf_c_fit_prev, nperm = 1)
crf_c_varimp_cond_prev <- varimp(crf_c_fit_prev, conditional = T, nperm = 1)

crf_c_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_c_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF prevalence ~ extraversion

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_e_fit_prev <- cforest(slope_prev ~ pers_e + women + academics +
                          cdu + afd + hospital_beds + gdp + manufact +
                          airport + age + popdens, 
                         df_ger_slope_prev[-1], 
                         controls = ctrls)

crf_e_varimp_prev <- varimp(crf_e_fit_prev, nperm = 1)
crf_e_varimp_cond_prev <- varimp(crf_e_fit_prev, conditional = T, nperm = 1)

crf_e_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_e_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF prevalence ~ agreeableness

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_a_fit_prev <- cforest(slope_prev ~ pers_a + women + academics +
                          cdu + afd + hospital_beds + gdp + manufact +
                          airport + age + popdens, 
                         df_ger_slope_prev[-1], 
                         controls = ctrls)

crf_a_varimp_prev <- varimp(crf_a_fit_prev, nperm = 1)
crf_a_varimp_cond_prev <- varimp(crf_a_fit_prev, conditional = T, nperm = 1)

crf_a_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_a_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF prevalence ~ neuroticism

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_n_fit_prev <- cforest(slope_prev ~ pers_n + women + academics +
                          cdu + afd + hospital_beds + gdp + manufact +
                          airport + age + popdens, 
                         df_ger_slope_prev[-1], 
                         controls = ctrls)

crf_n_varimp_prev <- varimp(crf_n_fit_prev, nperm = 1)
crf_n_varimp_cond_prev <- varimp(crf_n_fit_prev, conditional = T, nperm = 1)

crf_n_varimp_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_n_varimp_cond_prev %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF social distancing ~ openness

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_fit_socdist <- cforest(slope_socdist ~ ., 
                         df_ger_slope_socdist[-1], 
                         controls = ctrls)

crf_varimp_socdist <- varimp(crf_fit_socdist, nperm = 5)
crf_varimp_cond_socdist <- varimp(crf_fit_socdist, conditional = T, nperm = 5)

crf_varimp_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_varimp_cond_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF social distancing ~ conscientiousness

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_c_fit_socdist <- cforest(slope_socdist ~ pers_c + women + academics +
                          cdu + afd + hospital_beds + gdp + manufact +
                          airport + age + popdens, 
                         df_ger_slope_prev[-1], 
                         controls = ctrls)

crf_c_varimp_socdist <- varimp(crf_c_fit_socdist, nperm = 5)
crf_c_varimp_cond_socdist <- varimp(crf_c_fit_socdist, conditional = T, nperm = 5)

crf_c_varimp_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_c_varimp_cond_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF social distancing ~ extraversion

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_e_fit_socdist <- cforest(slope_socdist ~ pers_e + women + academics +
                          cdu + afd + hospital_beds + gdp + manufact +
                          airport + age + popdens, 
                         df_ger_slope_prev[-1], 
                         controls = ctrls)

crf_e_varimp_socdist <- varimp(crf_e_fit_socdist, nperm = 5)
crf_e_varimp_cond_socdist <- varimp(crf_e_fit_socdist, conditional = T, nperm = 5)

crf_e_varimp_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_e_varimp_cond_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF social distancing ~ conscientiousness

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_a_fit_socdist <- cforest(slope_socdist ~ pers_a + women + academics +
                          cdu + afd + hospital_beds + gdp + manufact +
                          airport + age + popdens, 
                         df_ger_slope_prev[-1], 
                         controls = ctrls)

crf_a_varimp_socdist <- varimp(crf_a_fit_socdist, nperm = 5)
crf_a_varimp_cond_socdist <- varimp(crf_a_fit_socdist, conditional = T, nperm = 5)

crf_a_varimp_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_a_varimp_cond_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

CRF social distancing ~ neuroticism

ctrls <- cforest_unbiased(ntree=500, mtry=5)

crf_n_fit_socdist <- cforest(slope_socdist ~ pers_n + women + academics +
                          cdu + afd + hospital_beds + gdp + manufact +
                          airport + age + popdens, 
                         df_ger_slope_prev[-1], 
                         controls = ctrls)

crf_n_varimp_socdist <- varimp(crf_n_fit_socdist, nperm = 5)
crf_n_varimp_cond_socdist <- varimp(crf_n_fit_socdist, conditional = T, nperm = 5)

crf_n_varimp_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))

crf_n_varimp_cond_socdist %>% as.data.frame() %>% rownames_to_column('variable') %>%
  ggplot(aes(x=variable, y=.)) +
  geom_bar(stat = 'identity') + 
  theme(axis.text.x = element_text(angle = 90))